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We employ the formalism of bond currents, expressed in terms of the nonequilibrium Green 
functions, to image the charge flow between two sites of the honeycomb lattice of graphene ribbons 
of few nanometers width. In sharp contrast to nonrelativistic electrons, current density profiles 
of quantum transport at energies close to the Dirac point in clean zigzag graphene nanoribbons 
(ZGNR) differs markedly from the profiles of charge density peaked at the edges due to zero-energy 
localized edge states. For transport through the lowest propagating mode induced by these edge 
states, edge vacancies do not affect current density peaked in the center of ZGNR. The long-range 
potential of a single impurity acts to reduce local current around it while concurrently increasing 
the current density along the zigzag edge, so that ZGNR conductance remains perfect G = 2e 2 /h. 
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Introduction. — The recent experimental discovery of 
a two-dimensional (2D) allotrope of carbon, termed 
graphene, has ushered unforeseen avenues to explore 
transport and interactions of low-dimensional electron 
system, build quantum-coherent carbon-based nanoelec- 
tronic devices, and probe high energy physics of "charged 
neutrinos" in table-top experiments Graphene rep- 
resents one-atom-thick layer of carbon atoms tightly 
packed into a honeycomb crystal lattice whose symme- 
tries impose linear energy-momentum dispersion of the 
low-energy quasiparticles [2]. Moreover, its bipartite 
structure introduces an internal pseudospin degree of 
freedom which connects electrons and holes through chi- 
rality (projection of pseudospin on the direction of mo- 
tion) of opposite signs, so that the effective mass equation 
turns into the Weyl equation for massless Dirac fermions 
(such as neutrinos) [2[. 

Relativistic energy spectrum, pseudospin, and zero 
gap with linearly vanishing density of states in the bulk 
graphene, have been probed in transport experiments un- 
veiling the 'chiral' quantum Hall effect, 'minimal conduc- 
tivity' at the charge neutrality (Dirac) point Ep =0, 
and weak-localization-type of quantum interference ef- 
fects [1]. The intriguing concept of chirality, whose con- 
servation would be responsible for the suppression ^ 
of backscattering from smooth (on the scale of the lat- 
tice constant) disorder and Klein tunneling Q through 
high and wide electrostatic potential barriers, has also 
led to a number of theoretical predictions [4| for esoteric 
micrometer-size graphene-based devices. 

On the other hand, recent experiments on graphene 
wires of nanoscale width have demonstrated the existence 
of a gap in their energy spectrum, Q which would allow 
GNR to replace semiconductor single-wall carbon nan- 
otubes while allowing for an easy integration into nano- 
electronic circuits via standard lithography end etching 
techniques. Direct STM imaging of the states localized 
at the edges of realistic GNR [6] , as well as possible chi- 
rality non-conserving scattering off the GNR edges, re- 



quires to examine the effects of edge-topology-dependent 
transverse subband structure [7, 8], edge states, impuri- 
ties, and potential barriers in tailoring quantum trans- 
port properties of GNR-based devices. 

The effect of zero-energy quantum states localized at 
the edges of ZGNR shown in Fig. [T] (which originate 
from the gauge field generated by lattice deformation 
and reflect the topological order in the bulk of bipartite 
honeycomb lattice [1Q[), as well as the energy gaps in 
armchair graphene nanoribbons (AGNR) controllable by 
their width, have been studied theoretically for more than 
a decade in equilibrium [111 Il2| and conduction proper- 
ties @, H, EH, 13, Gil However, very little is known 
about local features of transport through GNR. Fur- 
thermore, the application of recently advanced scanning 
probe techniques, developed to image local charge flow 
in quantum transport through 2D electron gases buried 
inside semiconductor heterostructures |17|, to graphene 



samples is eagerly awaited [lj . Many interesting findings 
are anticipated [l| when 2D electron states exposed on 
graphene surface are directly accessed by tunneling and 
local probes. Also, such transport experiments, going 
beyond traditional measurements of macroscopically av- 
eraged quantities, are becoming increasingly important 
for the development of nanolectronic devices — for exam- 
ple, recent imaging [18[ of the charge flow in conventional 
p — n junctions suggests that in structures shrunk below 
50 nm individual positions of scarce dopants will affect 
their function, thereby requiring to know precisely how 
charge carriers propagate on the nanoscale. 

Here we extend the bond current formalism for square 
lattices H, 0, EH to graphene honeycomb lattice, which 
allows us to predict spatial profiles of nonequilibrium 
charge and current densities. This imaging of charge flow 
provides direct insight into how massless Dirac fermions 
propagate between two neighboring lattice sites. The 
two-terminal device setup is shown in Fig. [T] where finite 
GNR sample is attached to two semi-infinite GNR leads. 
When the sample is clean, the whole structure represents 
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FIG. 1: (Color online) (a) Two terminal device, biased by the 
electrochemical potential difference \il — \ir — eV, consisting 
of a finite ZGNR sample connected to two ideal semi-infinite 
ZGNR leads. The width of 10-ZGNR is measured by the num- 
ber of zigzag chains N z = 10, and a is the lattice constant. 
The definition of bond current between two sites J mm / , shown 
in panel (a) by arrows connecting two supercells, and outflow- 
ing current J™* at site m is illustrated in panel (b). Panel 
(a) also plots the pattern of the local density of states at 
Ef = 0.0 It, which is dominated by the zero-energy localized 
edge states, as shown in panel (c) and in the corresponding 
total density of states in its inset. 



infinite ZGNR (illustrated by Fig.[Q) or AGNR, while the 
disordered sample is created by introducing vacancies at 
its edges or short-range or long-range impurity potential 
within its interior. The whole structure is described by 
the tight-binding Hamiltonian 

H — ^ ^ ^mm'^Cm' H~ ^ ^ ^m^m^m? (1) 
mm' m 

with single 7r-orbital per site. Here c m (c m ) creates 
(annihilates) an electron in the 7r-orbital at the site 
m = (m x ,m y ), and t mm / = t ~ 2.84 eV is a hopping 
parameter between nearest neighbor orbitals (which al- 
lows to reproduce ab initio [22[ computed structure of 
the conduction and valence bands in the vicinity of K 
and K' Dirac points located in two inequivalent corners 
of the hexagonal Brillouine zone where the bands touch 
conically). The impurity potential at site m is V m . 

Bond current formalism on graphene honeycomb 
lattice. — The central quantity of the steady-state lo- 
cal transport formalism on tight-binding lattices is the 
nonequilibrium lesser Green function G^ m (r = 0) = 
t(4^m> = U^dEG^E) function 0, where 
(...) is the nonequilibrium statistical average with re- 
spect to the density matrix at time t' = 0. It yields the 



magnitude of the bond current 

E F +eV/2 

Jmm' = J dE ^ t m i m G m m / {E) t m m i G m / m (-£/)] , 

Ep-eV/2 

(2) 

between the lattice sites m and its nearest neighbor site 
m', and the nonequilibrium charge density at site m 

E F +eV/2 
E F -eV/2 

These are the expectation values of the corresponding 
operators, J mm / = ^J mm ^ and n m = e^7V m ^, which 
satisfy the charge continuity equation on the lattice, 
edN m /dt + J2 m > (jmm' - Jm'm^j = 0, where m' are the 
three nearest neighbor sublattice B sites of site m be- 
longing to sublattice A, and vice versa. Thus, the bond 
current J mm / can be visualized as a bundle of flow lines 
bunched together along a link joining the two sites. 

The connection between the bond current vector J mm / , 
outflowing current from site m and total current 

1 = E m ' l J mm'| = E m ' J ™m> is illustrated by Fig.HKb). 

i i 

That is, when magnitudes of all vectors Jmm' (where 
length of the arrow is proportional to local current) on 
the bonds connecting supercells in Fig. [Tfa) are summed 
up to get /, then G = I /V gives the linear response con- 
ductance for small applied voltage bias V = O.OOlt/e. 
For zero-temperature quantum transport of electrons in- 
jected at the Fermi energy Ep = O.Oli there is only one 
open conducting channel, so that G = I /V = 2e 2 /h = 
Gq (Gq is the conductance quantum) for spatial distri- 
bution of local currents within ZGNR of Fig. [TJa). 

The matrix G < (E) contains information about the 
occupied states in the central region, and can be ob- 
tained by solving the Keldysh equation G < (E) — 
G r (E)"E < (E)G a (E). In the single-particle approxima- 
tion, where interactions are of the mean-eld type, this 
equation can be solved exactly by evaluating the re- 
tarded G r (E) = [E - H - U m - E£ - E^]- 1 and 
the advanced G a (E) = [G r (£')]"'" Green function ma- 
trices. In the absence of inelastic scattering, the re- 
tarded self-energies T, r L (E - eV/2) and T, r R (E + eV/2) 
introduced by the left and the right lead [24|, respec- 
tively, determine = -2i[Im5] L (£ - eV/2)f L (E - 
eV/2) + ImS^E + eV/2)f R (E + eV/2)} [f LyR are the 
Fermi distribution functions of the electrons injected 
from the macroscopic reservoirs through the leads and 
ImT, l^r = (T> r L R — T, a L R )/2i]. The gauge invariance of 
measurable quantities with respect to the shift of electric 
potential by a constant is satisfied on the proviso that 
J2 r R depend explicitly on the applied bias voltage V 
while G r (E) has to include the electric potential land- 
scape U m within the sample. 
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FIG. 2: (Color online) Spatial profiles of charge density n m 
and local current magnitude J™ 1 at each site of the cross 
section of ideal (disorder and defect free) infinite ZGNRs, 
AGNRs, and conventional quantum wires defined on square 
tight-binding lattices. The panels in the first row correspond 
to zero-temperature single channel quantum transport with 
the Landauer conductance G = Gq, while in the second row 
the Fermi energy Ef is tuned to allow for multichannel trans- 
port with G = 5Gq (Gq = 2e 2 /h). 



The key issue in the study of transport through nano- 
electronic devices attached to graphene leads is efficient 
algorithm to compute the retarded surface Green func- 
tion at the terminating edge of the semi-infinite GNR. 
For this purpose we employ the Ando algorithm [251 ], 
which constructs this Green function from the left and 
right transverse propagating exact Bloch eigenmodes. 
The original algorithm [25[ has to be generalized [26[ to 
handle non-invertible Hamiltonian Hi connecting super- 
cells consisting of atoms described by the Hamiltonian 
Ho (see Fig. [1] for illustration). 

Imaging charge flow in clean GNRs. — In ZGNR local- 
ized edge states appear at energies close to the Fermi 
level Ef — of undoped graphene. Thus, they mani- 
fest as edge peaks in the local density of states (DOS) 
p(E,m) = -lmG r mm (E)/7r in Fig. Q] at E = O.Olt, as 
well as a peak in the total DOS D(E) = Yj m P(E,m) 
at E = [where in the bulk graphene D(E = 0) = 0] 
shown in the same figure. They correspond to non- 
bonding molecular orbit als, where electrons are strongly 
localized near the zigzag edge composed on sublattice A 
sites (bottom) or sublattice B sites (top), and, therefore, 
cannot carry current. However, the overlap of two edge 
states from the top and bottom edge of a ZGNR yields 
bonding and anti-bonding states ensuring single conduct- 
ing channel (with partially flat energy-momentum dis- 
persion [7|, l8[) at energies arbitrarily close to the Dirac 
point. Therefore, ZGNR are metallic for all widths (as 
long as ferromagnetic ordering on zigzag edges is not 
taken into account [Hj]). This channel, together with 
2n right moving propagating modes crossing the Fermi 
energy for < \Ep\ < t yields odd- number conductance 
quantization G(E F ) = (2n + 1)Gq [H, 0, EH in clean 
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FIG. 3: (Color online) The linear-response (eV = O.OOlt) 
local current J^ 1 at each site of 10-ZGNR with single vacancy 
on the bottom zigzag edge. The current is proportional to 
the length of the arrow. In (a) the electrons are injected 
through a single conducting channel at Ef = O.Olt, while 
in (b) the Fermi energy of electrons injected through ZGNR 
semi-infinite leads is set to Ef = 0.8 at which there are five 
open channels. 



ZGNR. 

The Fermi energy of undoped GNR is at half-filling 
Ef = due to perfect electron- hole symmetry. In narrow 
ZGNR the gap A12 ~ 1/N Z between the second subband 
and E F = is so large (e.g., A 12 = OAt for N z = 10) that 
transport would remain within the single channel regime 
Ef < A12 for presently achievable densities of extra car- 
riers that can be injected into the ZGNR. Therefore, we 
focus on single-channel transport at Ef — O.Olt with 
maximum conductance G = Gq in the rest of the paper, 
while also using multichannel transport (Ef = O.St) with 
maximum conductance G = 5Gq for comparison. 

Figure [2] contrasts spatial profiles of and n m for 
single channel-transport through 10-ZGNR, 20-AGNR, 
and quantum wire modeled on a conventional square 
tight-binding lattice with 10 sites (hosting single s- 
orbitals) over its cross section. Following previous con- 
vention, the width of A/^-ZGNR is measured by the num- 
ber of zigzag chains N z , while 7V a -AGNR have N a dimer 
lines across the ribbon width. In Fig. [21 10-ZGNR corre- 
sponds to 20-AGNR and square lattice wires with N=10 
sites per cross section in the sense that all three quantum 
wires support maximum of 10 open conducting channels. 

The charge density (Fig. [2]) of low energy states is pro- 
portional to n(y) ex |</>a| 2 + |</>b| 2 , which in ZGNR (Fig.CQ) 
is peaked at its edges. On the other hand, the current 
density j x (y) v(&a x <$>) [a x is the Pauli matrix rep- 
resenting the x-component of the pseudospin operator 
acting on the AB space] of Dirac fermions, which are 
described in the low-energy (i.e., long wavelength) limit 
by continuum theory for the two-component wave func- 
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FIG. 4: The Landauer conductance of ZGNR with single 
impurity positioned in the center of the ribbon, to generate 
short-range d = 0.01a or long-range d = 1.5a Gaussian po- 
tential, as a function of the impurity potential strength. The 
electrons are injected from ZGNR lead through a single con- 
ducting channel (Ef — 0.01£) in the main panel and through 
five conducting channels (Ef — O.St) in the inset. 



tions & = (cj)* A (j)* B ) [defining relative contribution of the 
A and B sublattice in the make-up of quasiparticles] , 
reaches maximum in the center of the ribbon. This, to- 
gether with imaging of local charge flow J^* shown in 
Fig. 03 explains why isolated edge vacancies have very 
little effect on the Landauer conductance G « 0.98Gq 
of single- channel transport through ZGNR [l5[. When 
more channels are open for zero-temperature quantum 
transport, the vacancy affects not only the local current 
at the zigzag edge as in Fig. [31(a) , but also the magni- 
tude and the direction of within the bulk of ZGNR 
in Fig. [31(b), so that its conductance drops from G = 5Gq 
(when the vacancy is absent) to G « 4Gq. 

In contrast to ZGNR, spatial profiles of charge and 
current density in AGNR are highly inhomogeneous in 
both single and multichannel transport regimes. It is 
also instructive to compare n m vs. J™ 1 profiles in GNRs 
with those of non-relativistic electrons in quantum wires 
defined on the square lattice (third column in Fig. [2]). 
Their scalar wave function <j>(x,y) = </>trans(2/)e z/c:E yields 
the charge density n(y) cx \(j)(y)\ 2 which has the same 
profile as the corresponding current density j x (y) cx 
</)d x <f>* - <j>*d x <j> cx kn(y). 

Imaging charge flow in disordered GNRs. — In recent 
intense efforts to understand experimentally observed 
properties of the conductivity of bulk graphene (such as 
its linear dependence on carrier concentration, minimal 
value ~ e 2 /h at the Dirac point, and absence of weak 
antilocalization, expected due to chiral nature of elec- 
trons in graphene, or suppression of standard weak lo- 
calization corrections [l|), the analysis of scattering of 
massless Dirac fermions from different types of impu- 
rities has played an essential role [27j]. Since energy- 
momentum dispersion e(k x ) of transverse propagating 
modes in GNRs strongly depends on the confinement 
effects and topology of their edges, the investigation 
of the disorder effects in transport properties of meso- 
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FIG. 5: (Color online) The difference (disorder) - 
Jm* (clean) between the local currents at each site of a clean 
10-ZGNR and the same ZGNR with disorder introduced as a 
single impurity positioned on the sublattice A or sublattice B 
site in its center. The Gaussian impurity potential of strength 
Uo — 5.0t is short-range in (a) d — 0.01a and long-range in 
(b) d = 1.5a. The value of the corresponding conductances 
of single- channel (Ef = O.Olt) and five-channel (Ef — 0.8t) 
quantum transport is marked by arrows in Fig.2]and its inset, 
respectively. 



scopic graphene structures requires to reexamine condi- 
tions for the absence of backscattering due to conserva- 
tion of chirality [2]. For example, it has been argued 
recently M that large momentum difference between two 
valleys [8[ of ZGNRs at k x = 2ir/3a and k x = — 2ir/3a 
(which originate from the Dirac points K and K' of bulk 
graphene) prevents inter-valley scattering by long-range 
disorder, so that the special conducting channel gener- 
ated by the edges states remains a 'chiral' propagating 
mode in the sense that ZGNR has perfect conductance 
G(E F ) = 2e 2 /h for \E F \ < A 12 that does not decay 
as the length of the wire is increased [7] This special 
channel is comprised of states belonging to only one val- 
ley (unlike higher energy modes where both valleys con- 
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tribute), switching to the opposite one when changing the 
direction of propagation and allowing for valley filter de- 
vices [8]. The Landauer- formula-computed conductance 
does decay exponentially due to Anderson localization in 
the single or multichannel transport regimes when impu- 
rity potential is short-ranged [7]. 

To develop a real-space local picture of conduction 
through the special channel of ZGNR, we employ the 
Gaussian potential V m = U exp(— |m — m.o\/d 2 ) of range 
d generated by a single impurity centered at site mo, 
which belongs to either sublattice A or B. The poten- 
tial strength Uq is defined by the normalization condi- 
tion [3], J2 m U exp(— |m - m \/d 2 ) = Uq. The conduc- 
tance of ZGNR as a function of Uq for both short-range 
d — 0.05a (i.e., Anderson- model- type of disorder) and 
long-range d = 1.5a potential (which is short-ranged on 
the scale of the system size but varies smoothly on the 
atomic scale, corresponding to, e.g., screened charges in 
the substrate) is shown in Fig. [H in the single-channel 
(Ef = O.Olt) and multichannel (Ep = O.St) quantum- 
coherent transport regime. For long-range impurities, the 
conductance remains perfect G = 2e 2 /h up to Uq ~ 0.2t 
(and G ~ 1.99e 2 //i up to Uq ~ St), which is comparable 
to the energy difference between the transverse subbands. 

The corresponding images of local charge flow through 
ZGNR are shown in Fig. [5] by plotting the difference 
between local currents at each site (disorder) — 
J ^ (clean) , where Uq = in the clean case and Uq = bt 
in the disordered ZGNR. The arrows of local currents 
pointing to the right in Fig. [5] signify the reduction of 
current density due to the presence of impurity. While 
their sum in Fig. 03(a) indeed explain how current den- 
sity gets reduced around the short-range impurity, in the 
case of long-range impurity the same reduction of current 
density around the impurity is compensated by the in- 
crease of current density along the zigzag edge displayed 
in Fig. 03(b). Moreover, the edge being exploited by mass- 
less Dirac fermions to resist current degradation consists 
of the same sublattice sites as is the site on which the 
impurity is located. 

Conclusions. — We have shown how to adapt the bond 
current formalism to graphene honeycomb lattice which 
makes it possible to predict spatial profiles of local cur- 
rents of massless Dirac fermions between two sites of the 
lattice. The profiles we obtain for graphene nanoribbons 
suggest several interesting experiments that are within 
the reach of present local probe-based direct charge imag- 
ing techniques (i) in ZGNR with localized edge 
states the large part of current flows through the cen- 
ter of the ribbon which makes single-channel quantum 
transport largely insensitive to edge vacancies; (ii) while 
both short-range and long-range impurities reduce cur- 
rent density in the region of their influence, in the single- 
channel transport through the lowest transverse prop- 
agating mode generated by the edge states of ZGNR 
this reduction in the case of long-range impurities can 



be compensated by the increase of current density along 
the zigzag edge ensuring perfectly conducting channel 
G = 2e 2 /h even in the presence of disorder. 

We thank E. Andrei, A. H. MacDonald, and P. Kim 
for valuable discussions. 
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